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Abstract 

A numerical algorithm for solving mantle convection problems with strongly variable 
viscosity is presented. Equations for conservation of mass and momentum for highly 
viscous and incompressible fluids are solved iteratively by a multigrid method in 
combination with pseudo-compressibility and local time stepping techniques. This 
algorithm is suitable for large-scale three-dimensional numerical simulations, be- 
cause (i) memory storage for any additional matrix is not required and (ii) vector- 
ization and parallelization are straightforward. The present algorithm has been in- 
corporated into a mantle convection simulation program based on the finite- volume 
discretization in a three-dimensional rectangular domain. Benchmark comparisons 
with previous two- and three-dimensional calculations including the temperature- 
and/or depth-dependent viscosity revealed that accurate results are successfully re- 
produced even for the cases with viscosity variations of several orders of magnitude. 
The robustness of the numerical method against viscosity variation can be signifi- 
cantly improved by increasing the pre- and post-smoothing calculations during the 
multigrid operations, and the convergence can be achieved for the global viscosity 
variations up to 10^'^. 

Key words: mantle convection, variable viscosity, pseudo-compressibility, 
multigrid method, local time-stepping 



1 Introduction 



The Earth's mantle is the spherical shell composed of silicate rocks, and it 
ranges from approximately 5-50km to 2900km depth. Although the mantle 
behaves like an elastic solid on short time scales, it acts like a highly viscous 
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fluid on long time scales. The mantle also acts as a heat engine, and it is in a 
convective motion in order to transport the heat from the hot interior to the 
cool surface [1,2]. The mantle convection is observed as the motion of tectonic 
plates at the Earth's surface. The motion of surface plates, in turn, drives 
seismicity, volcanism and mountain building at the plate margins. Thus, the 
mantle convection is the origin of the geological and geophysical phenomena 
observed at the Earth's surface. A major tool for understanding the mantle 
convection is numerical analysis. It has been playing an important role in the 
study of mantle convection, since a numerical simulation of mantle convection 
first arose [3,4]. 

Mantle convection requires different numerical techniques from those for or- 
dinary fluids such as water because of its rheological properties. The viscosity 
of mantle materials is estimated as high as 10^^ Pa s [1.5]. Since the man- 
tle materials is highly viscous, both the nonlinear and time-derivative terms 
of velocity can be ignored in the equation of motion. This implies that the 
flow in the mantle is described by a steady-state Stokes flow balancing among 
the buoyancy force, pressure gradient and viscous resistance. Taken together 
with the assumption of incompressibility, one needs to solve elliptic differen- 
tial equations for velocity and pressure at every timestep. In addition, the 
viscosity of mantle material varies by several orders of magnitude depending 
on temperature, pressure, and stress [6,7]. The strong variation in viscosity 
makes numerical techniques for ordinary isoviscous fluids, such as the spectral 
method [8,9], unfit for the numerical modeling of mantle convection. In order 
to get deep insights into the mantle convection, it is very important to de- 
velop efficient numerical techniques that can deal with the steady-state flow 
of highly viscous and incompressible fluids with a strongly variable viscosity. 

The efficiency of numerical simulations of mantle convection strongly relies on 
numerical methods used for solving elliptic differential equations. One of the 
most efficient methods is the multigrid iteration [10]. The multigrid concept 
has been successfully applied to a wide range of problems, including calcu- 
lations of incompressible fluid flow [11,12,13]. During the last two decades, 
various numerical models of mantle convection have been developed where the 
multigrid method is utihzed. There are two strategies to apply the multigrid 
method to this problem, depending on how the steady-state Stokes equations 
are solved. 

The first strategy solves the Stokes equations by splitting into the separate 
equations for velocities and pressure. The discretized equations for velocity 
components (or their proxy) are solved by the multigrid method, while the 
pressure is ehminated or solved separately. Parmentier et al. [14] developed 
convection models of isoviscous fiuid in three-dimensional Cartesian geometry. 
By using a streamfunction formulation, the Stokes equations are reduced to a 
pair of Poisson equations which are solved by multigrid iterations. Baumgard- 
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ner and his colleagues developed convection models for fluids with constant 
[15] and moderately variable viscosity [16,17] in a three-dimensional spherical 
geometry. They solved the elliptic equations for velocity components using a 
multigrid method, while the pressure fields are prescribed by the equation of 
state. Moresi and his colleagues [18,19] developed a model for convection prob- 
lems with strongly variable viscosity in two- and three-dimensional Cartesian 
geometry. The Stokes equations are solved separately for velocity and pressure 
by so-called Uzawa iterative scheme. The iteration for velocity is carried out 
by a multigrid method, while a conjugate gradient scheme is used for pressure 
iteration. 

The second strategy, on the other hand, solves the Stokes equations for veloc- 
ity and pressure as a whole by the multigrid technique. The key issue of this 
strategy is a choice of the smoothing algorithm which reduces the errors of 
solution on a particular grid. Several methods for solving incompressible fluid 
flows have been utilized as a smoothing algorithm. Trompert and Hansen 
[20,21] and Albers [22] developed numerical methods for convection problems 
with variable viscosity in three-dimensional Cartesian geometry. The Stokes 
equations are solved by a multigrid method where the SIMPLER algorithm 
[23] is employed as a smoothing operation. Auth and Harder [24] used the 
symmetric coupled Gauss-Seidel (SCGS) method [25] as a smoothing opera- 
tion in the multigrid method for solving the Stokes equation. Tackley [26,27] 
employed a similar smoothing method where the velocity and pressure fields 
are updated in a somewhat coupled manner, and his method has been success- 
fully applied to a variety of convection problems [27,28], including calculations 
with extremely variable viscosity [29,30,31]. 

In this paper, we present a solution algorithm for the three-dimensional man- 
tle convection, following the second strategy where the Stokes equations are 
solved as a whole by the multigrid method. In Section 2 we introduce our so- 
lution algorithm based on the pseudo-compressibility method [32]. A notable 
feature of our algorithm is its simplicity and intuitiveness. This feature makes 
the algorithm easily fit to a smoothing operation of the multigrid iterations. In 
addition, taken together with the local time-stepping method, a strong vari- 
ation in viscosity can be handled without a severe increase in computational 
costs. In Section 3, we develop a numerical model of mantle convection that 
can deal with a strongly variable viscosity based on the algorithm presented 
in Section 2 together with the multigrid technique. In order to demonstrate 
the validity and efficiency of our algorithm, we carry out several calculations 
for the mantle convection with a strongly variable viscosity. In this paper, an 
application of our algorithm is the thermal convection in a three-dimensional 
rectangular domain. It can be easily applied to the numerical models for more 
realistic problems of mantle convection, such as spherical shell geometry and 
thermo-chemical convection and so on. 
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2 Iterative Algorithm for mantle convection 



We consider thermal convection of a highly viscous and incompressible fluid 
with a strongly variable Newtonian viscosity in a three-dimensional Cartesian 
geometry {xi—x, X2—y, x^—z). The nondimensional forms of the fundamental 
equations are [33,34, for example]; 



d 



dxj L 



^ dxj dxi 



dp 
dxi 
dvk 
dxk 

or 



or 

dt ' dx. 



+ ^ = RaTSi3 = 1,2,3), 

(A; = 1,2,3), 



(1) 
(2) 
(3) 



where Vi {i — 1,2, 3) are fluid velocities in i-th direction, p pressure, r) viscosity, 
T temperature, Ra the Rayleigh number, and q is the internal heating rate. 
We assumed that x^-axis is the vertical axis pointing upward. Boussinesq 
approximation is employed in the energy equation (3) and, hence, the effects 
of adiabatic and viscous heating are ignored. 



2. 1 Overview of the iterative algorithm 



In the following, we consider the solution algorithm for the conservation equa- 
tions of momentum (1) and mass (2). In order to simplify the solution algo- 
rithm, the energy equation (3) is solved separately from Eqs. (1) and (2): the 
temperature T at new time is calculated using the velocity at old time in the 
advection term. This algorithm solves Eqs. (1) and (2) for the velocity and 
pressure at new time using the new temperature. In addition, the viscosity r] 
at new time is already known by some means and kept unchanged while the 
velocity and pressure are computed. Because of these assumptions, Eqs. (1) 
and (2) are taken to be linear with respect to velocity and pressure. 

The discretized equations of (1) and (2) can be written in a matrix form of 
Ax = b, (4) 



where 



4 





Cn 


Ci2 








Vl 











^22 




G2 




V2 







A — 










, X = 




h — 








^32 




G3 




V3 




RaO 




.^1 











P . 








(5) 



In (5), (i = 1, 2, 3) are the vectors of unknown fluid velocities in i-direction, 

p the vector of unknown pressures, 6 the vector of known temperatures, Cij 
{i.j = 1. 2, 3) the matrices representing the spatial discretization of velocities 
for calculating viscous stress, Gj (i = 1, 2, 3) the matrices representing the 
discretization of the pressure gradient, {i — 1,2,3) the matrices represent- 
ing the discretization of the divergence of velocity, and is the zero matrix 
or zero vector. One of the major difficulties in solving (4) is that a zero block 
appears in the main diagonal in A. This implies that basic iterative methods 
for solving linear simultaneous equations (such as Jacobi and Gauss-Seidel 
methods) which make use of the inverse of the main diagonal fail to solve (4) . 

The algorithm for solving (4) that we propose in this paper is another kind of 
iterative method: Let x"- be an approximate at n-th iteration (n = 0, 1, 2, • • •). 
We then calculate the approximate at the iteration (n + 1) by, 

a;"+i = a;" + T(6- Aa;'*). (6) 



Here a regular matrix T is introduced so as to control the convergence rate 
of (6), and is chosen to be a diagonal matrix whose nonzero elements are 
positive. This iterative procedure is used together with the multigrid method, 
and is repeated until a;" converges to a steady solution. It is obvious that 
the steady solution of (6) satisfies (4), and that it does not depend on T. In 
addition, the presence of zero block in the main diagonal of A does not spoil 
the convergence of this iterative procedure. 

This iterative algorithm is suitable for large-scale numerical simulations on 
massively parallel computers and, in particular, massively vector-parallel su- 
percomputers because (i) the procedure (6) consists of multiplication of matri- 
ces and vectors and summation of vectors and, hence, can be easily vectorized 
and parallelized, and (ii) large memory storage for the matrix T is not re- 
quired, since T is a diagonal matrix. 

Another benefit of the iterative procedure (6) is that it is easily implemented 
into the multigrid method. This is because the procedure (6) is quite similar 
to that representing the basic iterative methods which are commonly used as 
a smoothing operator of the multigrid method. Suppose an iterative procedure 
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for (4) by Jacobi method; 

+ Ace") (n = 0, 1, 2 • • •), 



(7) 



where a matrix D is the main diagonal of A. The two procedures (6) and 
(7) are the same except that the diagonal matrix T is used in (6) in place of 
in (7). Therefore, the iterative procedure (6) can be used as a smoothing 
operator of the multigrid method in a similar manner to Jacobi method. 

In the following subsections, we will introduce the ideas which led us to the 
iterative procedure (6). We will also discuss an appropriate choice of the matrix 
T, which is a key parameter in this procedure. 



2.2 Ingredient 1: pseudo-compressibility method for highly-viscous, incom- 
pressible fluids 



When the matrix T is diagonal, Eq. (6) can be explicitly written as {i,j, k 
1,2,3), 



+ 



dxi dx-i 



\ dxi dxj J J 



Ra T5, 



i3 ( , 



^dxk 

where r^. (i = 1, 2, 3) and Tp are the diagonal elements appearing in T. 



(8) 
(9) 



In order to discuss the physical meaning of these equations more clearly, we 
introduce the following "auxiliary" set of equations; 



' dt 
dt 



dp 



d 



dxi dxj 
dvk 



dvi dv.: 



dxi dx4 



+ RaT5i3, 



dxk 



(10) 
(11) 



Here t is analogous to time, M^^ {i = 1,2,3) and K are positive constants 
analogous to density and compressibility, respectively. These equations come 
from the modification of the "auxiliary" set of equations used in the pseudo- 
compressibihty method [32] , which is used to solve steady-state flow of incom- 
pressible fluid with high Reynolds number. The difference between Eqs. (10) 
and (11) and those employed for high Reynolds- number flow is that the non- 
linear term of velocity has been eliminated in (10), because the viscosity of 



6 



mantle materials is significantly high and, in other words, the Reynolds num- 
ber is significantly small. 

By discretizing (10) and (11) in the direction of t by a first-order explicit 
scheme, we get 



n+l _ n 



K 



At 



dp"" d 
+ 



dxi dxj 



'dvl dvf 



At 



dxk 



(12) 
(13) 



We notice that Eqs. (12) and (13) and Eqs. (8) and (9) are identical if we 
choose Ty. and Tp as. 



r^. = Ai/M^., Tp = Ai/K. 



(14) 



In other words, r^. and Tp are "effective" timesteps for the evolution of Vi and 
p, respectively. 

The iterative procedure given by (12) and (13) converges to an incompressible 
fiow field regardless of the choice of My. and K. To show this more clearly, 
we here assume that M^-, K and 7] are constant. From Eqs. (10) and (11) we 
obtain the pseudo-temporal evolution of V • as, 

£l(V . v) = -^VHV ■ v) + [V^ (V • v)] . (15) 
dp ^ ' KM ^ ' Mdt^ ^ ^ ' 

(Here we denote My^ = M .) The first term of the right-hand side represents the 
effect of propagation of "pseudo-sound wave" whose velocity is 1/ VaM, and 
the second term represents the effect of diffusion due to viscosity. This equation 
indicates that the pseudo-temporal evolution of V • v is characterized by a 
decaying oscillation. Therefore, V • v approaches to an asymptotic value (= 0) 
as t increases to infinity. Consequently, the procedure of (12) and (13) leads 
the steady velocity field with V • t; = 0, as long as the numerical integration 
scheme is stable. 

We also note that the iterative procedure (8) and (9) should be used in combi- 
nation with multigrid method [10]. Recall that we are trying to find a steady- 
state solution of evolution equations (10) and (11) through a repetition of 
temporal integration, which inevitably requires a large number of iterations. 
Moreover, the elliptic nature of (10) for velocities Vi results in a slow re- 
duction in their errors, particularly in their smooth components with large 
spatial wavelengths. By incorporating into the multigrid procedure, we are 
able to obtain a sufficiently fast convergence. In fact, the multigrid technique 
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is adopted in most of recent numerical analysis using pseudo-compressibility 
method [35,36, for example]. 



2.3 Ingredient 2: local time stepping method for strongly variable viscosity 



In most of previous numerical analysis using pseudo-compressibility method 
[35,36,37,38, for example], the viscosity of fluids has been assumed to be con- 
stant. In this subsection, we apply this method for the cases with a strong 
variation in viscosity. 

The convergence rate of the original pseudo-compressibility method deterio- 
rates in proportion to the viscosity variation in the entire domain, when the 
spatial variation in viscosity is introduced. This is because the viscosity is a 
diffusion coefficient for Vi in Eq. (10). The spatial variation in rj results in the 
spatial variation in the local convergence rate of Vi. The values of Vi in the re- 
gion with smaller t] reach to an appropriate solution only very slowly, whereas 
the values in the region with larger 77 converge quickly. This suggests that the 
convergence should be accelerated particularly in the region with small rj. 

Here we try to accelerate the convergence for the case with a spatial variation 
in 77 by controlling the stepping of the pseudo-time t in accordance with r). 
This idea is known as the local time stepping method, which is used to obtain 
steady-state solution of evolutionary equations. This method allows updating 
each variable using a timestep which is based on the local numerical stability 
criterion. We thus spatially vary the effective timesteps t„. and Tp in (8) and 
(9) in accordance with 77. 

The spatial variation of r^,. and Tp can be estimated from the local numerical 
stability criterion of (15). In the following derivation, we assume for simplicity 
that Eq. (15) holds even in the case of variable viscosity. Prom the criterion 
of the diffusion of V ■ v (the second term in the right-hand side of (15)) we 
require that the nondimensional parameter a defined by 

_ 277/M 

a = ^ ^ , ^ ~ (16) 



is sufficiently small. Here A is the mesh size. Similarly, from the criterion of 
the pseudo-sound propagation we require that the nondimensional parameter 
/3 defined by 



(17) 

^ A/At ^ ' 
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is sufficiently small. By substituting (16) and (17) into (14) we get 

At aA^ At 2/52 

Here we assumed My. = M. 

When Ty^ and Tp are defined by (18), the spatial variations of the diffusion of 
velocity components and the pseudo-sound propagation become more modest 
than that of viscosity. Equation (18) implies that M (= M„. in (10)) and K 
are taken to be proportional to rj and r]^^, respectively. By assuming that 
ri/M and Krj are constant in (10) and (11), we obtain the pseudo-temporal 
evolution of V ■ v in the presence of spatial variation in 7] as. 



di 



72 (V-^^) 



1 



+ 



KM 

d_ 

KM dxi 



V^(V • v) + 



dxi 



2r] d 
Mdt 

{V-v) 



V^(V • v) 



+ 



d 



Mdxi 



dihirf) d 



dxj 



' dvj dvi 
dxi dxj 



(19) 



The first and second terms in the right-hand side of (19) are equivalent with 
those in the right-hand side of (15). Note that the effects of spatial variation in 
T) appear in the third and fourth terms in the form of In 77. Thus we expect that 
V • V reduces to zero in the entire computational domain at a rate comparable 
to that for the case with constant 77. 

Another constraint for t^^ and Tp can be obtained directly from (6). Let = 
x'" — A~^b and r" = b — Ax^ to be the error and the residual at n-th iteration, 
respectively. Prom (6), we get the evolution equations for e" and r" as. 



-,n+l 



[I - ATy 



(20) 



where / is the identity matrix. These equations indicate that the convergence 
rate of (6) becomes optimal when T = Topt = A~^. (In this sense, the ma- 
trix T can be regarded as a "preconditioner" for Eq. (4).) In the present 
algorithm using local time stepping, we approximate T ~ Topt by a diagonal 
matrix. In addition. Equation (20) implies that the necessary condition for 
the convergence is that the spectral radius of the matrix I — TA (or, identi- 
cally, I — AT) is smaller than one. Indeed, this condition is identical to the 
Courant-Friedrichs-Lewy (CFL) condition of Eqs. (12) and (13), i.e., a < etc 
and (3 < (3c where ac and Pc are threshold values. 
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3 Results 



The algorithm described in the previous section has been incorporated into a 
mantle convection simulation program in a three-dimensional rectangular do- 
main. We carried out several calculations for thermal convection with strongly 
variable viscosity, in order to demonstrate the validity and efficiency of the 
iterative algorithm. 



3.1 Details of numerical method 

The basic equations are discretized by the finite volume method. A staggered 
grid is used; temperature and pressure are located at the center of the grid 
cells, while the velocity components are at the center of the cell faces normal 
to the direction of the velocity components. A uniform mesh is employed. The 
basic equations are nondimensionalized with a length scale of h (height of the 
box), time scale oih? / n (where n is thermal diffusivity), temperature scale of 
AT (temperature drop across the box), and pressure scale of rjon/h'^ (where 
rjo is reference viscosity). In the following subsections, we present two kinds 
of calculations. First is to obtain the instantaneous flow flclds for prescribed 
distributions of buoyancy and viscosity, and second is to obtain the temporal 
developments of thermal convection. 

When the purpose is to obtain the flow flelds, we only solve the conservation 
equations for momentum (1) and mass (2). These equations are solved by the 
multigrid method based on correction-storage algorithm, because of the linear 
nature of (4). The smoothing operation at each grid level is carried out by (6). 
The diagonal elements of T for unknown velocities (r^.) were chosen to be 0.5 
times the inverse of the corresponding diagonal elements of A, while those for 
unknown pressures (r^) are to be 0.25 times the viscosity at the corresponding 
cell centers. The values of r^,- and Tp dcflncd thus are approximately equal to 

those with a = /3 = -g in (18). (We observed that the iteration by (6) diverged 
when Ty. or Tp is larger than the above values.) A linear interpolation is used in 
both flne-to-coarse (restriction) and coarse-to-flne (prolongation) operations. 
The values of viscosity at cell centers on coarser grids is calculated by a linear 
interpolation from those on the flner grids by one grid level. On the other 
hand, the values of viscosity at the midpoints of cell edges on each grid level, 
which are necessary to calculate viscous shear stress, are calculated by an 
interpolation proposed by Ogawa et al. [34] from the values at cell centers on 
the same grid level. Since the aim of this paper is to demonstrate the efficiency 
and robustness of the smoothing algorithm (6), only the V-cycle is used in 
the present calculations, although more complicated multigrid iteration cycles 
(such as W-, F-cycles) are expected to be more robust [24,22, for example]. 
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When the purpose is to obtain a time-dependent or steady-state convection, we 
solve the energy equation (3) in addition to Eqs. (1) and (2). In the calculations 
presented in this paper, the internal heating rate q is taken to be zero. The 
energy equation is discretized by a first-order Euler method in time. An upwind 
scheme, called power-law scheme [23] , is used to evaluate the contributions of 
heat transport by advection and conduction. The discretized equation for T is 
solved by a fully implicit scheme when we seek for a steady-state convection, 
or by a fully explicit scheme when we seek for a time-dependent convection. 

3.2 Benchmark comparison 

In order to test the validity of our numerical code, we compare our results 
with those of earlier studies [39,40,34]. The temperature T is fixed at and 
1 at the top {z — 1) and bottom {z — 0) boundaries, respectively, while the 
vertical side walls arc adiabatic. In most cases, the viscosity 77 is assumed to 
be dependent on temperature T and depth [1 — z) as, 

T] = 77texp[-£;TT + Ep{l - z)], (21) 

where r]t is the viscosity at the top surface {z — 1, T — Qi), and Et and Ep 
are constants describing the temperature- and depth-dependence of viscosity, 
respectively. 

3.2.1 Benchmark for two-dimensional convection 

First we present the benchmark comparison with Blankcnbach et al. [39] for 
the cases of steady-state convection in two-dimensional rectangular domains 
of aspect ratio (width/depth) A. We carried out calculations of five cases hsted 
in Table 1. Cases la to Ic are the convection of isoviscous fluid in a square 
box (A = 1), Case 2a is the convection of fluids with temperature-dependent 
viscosity in a square box, and Case 2b is the convection with temperature- 
and depth-dependent viscosity in a box of A = 2.5. The impermeable and 
shear-stress-free conditions are adopted along the all boundaries. The initial 
conditions are chosen to make single convective cell with an ascending flow 
along the left-side wall (x = 0). We carried out these calculations by varying 
the number of mesh divisions in x- and ^-directions, while the mesh division 
in y-direction was kept to 1, in order to obtain two-dimensional flow patterns 
in the x-z plane. 

We summarized in Table 1 the parameter values and mesh divisions employed 
in the present calculations, and the obtained values of (i) Nusselt number A^-u, 
(ii) root-mean-square velocity Kms, and (iii) vertical temperature gradient qi 
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Table 1 

Benchmark comparison with Blankenbach et al. [39] . See the text as for the meaning 



of the symbols. 



mesh 


16x16 


32x32 


64x64 


128x128 


benchmark standard 


(a) Case la (constant 


viscosity, Ra 


= lo^ A = 


1) 




Nu 


4.726885 


4.840006 


4.872745 


4.881429 


4.884409±0.000010 




41.755217 


42.560118 


42.785441 


42.844667 


42.864947±0.000020 


qi (= 93) 


7.695426 


7.966605 


8.036044 


8.053536 


8.059384±0.000003 


92 (= 94) 


0.887088 


0.664927 


0.608154 


0.593695 


0.588810±0.000003 


(b) Case 


lb (constant viscosity, Ra 


= lo^ A = 


1) 




Nu 


9.235550 


10.082608 


10.394786 


10.495490 


10.534095±0.000010 




181.18346 


189.19219 


191.98526 


192.87381 


193.21454±0.00010 


91 (= 93) 


14.424512 


17.764322 


18.731196 


18.990358 


19.079440±0.000040 


92 (= 94) 


3.441283 


1.425553 


0.898253 


0.767465 


0.722751±0.000020 


(c) Case 


Ic (constant 


viscosity, Ra 


= 10^ A = 


1) 




Nu 


13.505592 


18.963003 


20.884342 


21.604074 


21.972465±0.000020 


^ms 


675.69469 


780.80943 


814.68028 


827.43096 


833.98977±0.00020 


91 (= 93) 


16.63581 


31.21463 


41.44228 


44.72368 


45.96425±0.00030 


92 (= 94) 


9.836338 


6.941879 


2.477731 


1.255997 


0.877170±0.000010 


(d) Case 2a {Rat = 10^, Et = ln(103), Ep = 0, X = l) 




Nu 


10.3634 


10.4780 


10.1666 


10.0862 


10.0660±0.00020 




372.2759 


457.4738 


475.5292 


478.9647 


480.4334±0.1000 


91 


17.42173 


20.41169 


18.33101 


17.72910 


17.53136±0.00400 


92 


2.25932 


1.16135 


1.03954 


1.01498 


1.00851±0.00020 


93 


14.1505 


20.2221 


25.3614 


26.6611 


26.8085±0.0100 


94 


7.036967 


4.248769 


1.557456 


0.743177 


0.497380±0.000100 


mesh 


40x16 


80x32 


160x64 


320x128 


benchmark standard 


(e) Case 2b {Rat = 10^, Et = ln(16384), Ep = 


ln(64), A = 


2.5) 


Nu 


7.4020 


7.2761 


7.0450 


6.9605 


6.9299±0.0005 




201.478 


171.085 


175.343 


172.806 


171.755±0.020 


91 


18.4789 


30.6337 


20.9644 


19.0478 


18.4842±0.0100 


92 


0.37860 


0.17729 


0.18048 


0.17873 


0.17742±0.00003 


93 


14.5474 


14.7261 


14.7092 


14.3514 


14.1682±0.0050 


94 


3.51882 


1.28403 


0.79348 


0.66355 


0.61770±0.00005 
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Figure 1 



Fig. 1. Isotiiermal surfaces obtained for the benchmark calculations of stationary 
convections in Busse et al. [40]. (a) Case la is for constant viscosity, while (b) Case 
2 is for modestly temperature-dependent viscosity whose viscosity contrast is 20. 
The calculations were carried out with (a) 64 x 32 x 64 and (b) 64 x 64 x 64 mesh 
divisions. 

to ^4 at {x, z) = (0, 1), (A, 1), (A, 0), and (0, 0), respectively. We also show in 
Table 1 the results from the benchmark standard [39] . Table 1 shows that the 
agreement between the results of our calculations and benchmark standards 
is satisfactory for all the values. We conclude that the present numerical code 
accurately handles two-dimensional convection of fluids with both constant 
and variable viscosity. 



3.2.2 Benchmark for three-dimensional convection with modest viscosity vari- 
ation 

Second, we compare our results with those of Busse et al. [40] for the cases of 
steady-state convection with both constant (Case la) and modestly temperature- 
dependent viscosity (Case 2) in a three-dimensional rectangular box of axbx 1. 
Figure 1 shows the convective flow patterns for both cases. Case la is a bi- 
modal convection of isoviscous fluid with Ra — 3x 10^ in a box of a = 1.0079 
and b — 0.6283, while Case 2 is a square-cell convection in a cube {a — b — 1) 
with modestly temperature-dependent viscosity. In Case 2, the viscosity rj 
depends on temperature T as. 



ri{T) = exp 



Q Q 
TTT7 ~ 0.5 + 



(22) 



Here, the viscosity contrast is set to be 20, and the Rayleigh number deflned 
by the viscosity for T = 0.5 is 2 x 10^. In these cases the top and bottom 
surfaces are assumed to be rigid {vi = V2 = = 0), while the vertical side 
walls are planes of mirror symmetry. The initial conditions are chosen to make 
single ascending and descending flow at {x,y) = (0,0) and (a, 6), respectively. 

We summarize the result of benchmark comparison for three-dimensional con- 
vection in Table 2. We calculated the values of (i) Nu and Kmsj (ii) vertical 
velocity w and temperature T at specifled points at mid-depth of the con- 
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Table 2 

Benchmark comparison with Busse et al. [40]. See the text as for the meaning of 
the symbols. 

(a) Case la (constant viscosity, Ra = 3 x W^, a = 1.0079, 6 = 0.6283) 



mesh 


16x8x16 


32x16x32 


64x32x64 


128x64x128 


benchmark standard 


Nu 


3.7373 


3.6019 


3.5549 


3.5419 


3.5374±0.0005 




42.047 


41.382 


41.104 


41.026 


40.999±0.004 


u;(0, 0,0.5) 


122.919 


120.451 


117.880 


116.964 


116.625±0.030 


w{0, b, 0.5) 


20.044 


33.693 


38.566 


39.994 


40.500±0.030 


T(0, 0,0.5) 


0.79161 


0.80376 


0.80262 


0.80169 


0.80130±0.00005 


r(0,6,0.5) 


0.57314 


0.60324 


0.61434 


0.61761 


0.61876±0.00005 


Q(0,0) 


7.68656 


7.08517 


6.82042 


6.74082 


6.7127±0.0500 


Q{a,0) 


1.85083 


1.60300 


1.53359 


1.51460 


1.5080±0.0500 


Q{0,b) 


2.75670 


2.96542 


3.10916 


3.15671 


3.1740±0.0500 


Q{a,b) 


0.74943 


0.70104 


0.70778 


0.71225 


0.7140±0.0500 


(b) Case 2 (modestly temperature-dependent viscosity, a = b = 1) 




mesh 


16x16x16 


32x32x32 


64x64x64 


128x128x128 


benchmark standard 


Nu 


3.0870 


3.0503 


3.0419 


3.0399 


3.0393±0.0050 




35.350 


35.161 


35.130 


35.124 


35.13±0.05 


u;(0, 0,0.5) 


154.076 


162.020 


164.814 


165.631 


165.9±1.0 


w{0, b, 0.5) 


-25.662 


-26.446 


-26.650 


-26.703 


-26.72±0.1 


w{a, b, 0.5) 


-59.973 


-58.731 


-58.360 


-58.260 


-58.23±0.1 


T(0, 0,0.5) 


0.88286 


0.89803 


0.90323 


0.90474 


0.90529±0.0010 


r(0,6,0.5) 


0.51736 


0.50170 


0.49724 


0.49606 


0.49565±0.0010 


T(a, 6, 0.5) 


0.25785 


0.24417 


0.24052 


0.23957 


0.23925±0.0010 


g(o,o) 


6.09608 


5.90091 


5.85018 


5.83783 


5.834±0.015 


Q(a,0) 


1.80145 


1.73607 


1.71942 


1.71508 


1.714±0.015 


Q{a,b) 


0.78762 


0.77281 


0.76952 


0.76869 


0.768±0.015 



vecting vessel {z = 0.5), and (iii) vertical temperature gradient Q at specified 
points at the top surface {z — 0). We also show in Table 2 the results from 
the benchmark standard [40]. The agreement between the obtained results 

and the benchmark standard is satisfactory for all of the values. In particular, 
we obtained a good agreement for Case 2. We thus conclude that the present 
numerical code accurately handles three-dimensional convection of fluids with 
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Table 3 

Comparison of the flow patterns, Nussclt numbers Nu and the minimum mesh sizes 
in z-direction dziain between the present method ("KKS") and Ogawa et al. [34] 
("OSZ"). Note that in OSZ [34] the mesh spacing was taken to be non-uniform 
with finer resolution near the boundaries so as to improve the accuracy in Nu. 



Case 


Rat 




pattern 


code 




KKb 




OSZ 










mesh 


32x16x32 


64x32x64 


128x64x128 


24x14x22 










f]z ■ 


0.03125 


0.015625 


0.0078125 


3.2 X 10-2 


-I 
1 


11) 


i 


WL-3U 


l\u 


y.oi6 


i0.07i 


10.101 


n TO 
9.72 


4 


10 


1 a2 
10 


TTTT OTA 

WL-2D 


Nu 


4.060 


4.073 


4.077 


4.13 


10 


10^ 


10^ 


WL-3D 


Nu 


4.961 


4.932 


4.921 


4.96 


Case 


Rat 




pattern 


code 




KKS 




OSZ 










mesh 


ozx ibXoz 


d4XozXd4 


1ZoXd4X 12o 


/i 1 O on* 

44X loXoU 












0.03125 


0.015625 


n nnToi occ 
0.0078125 


1.7 X 10 


16 


10^ 


6.2 X 10 


TS7T • ) T ~\ 

WL-3D 


Nu 


5.284 


5.271 


5.271 


5.37 


Case 


Rat 


rr, 


pattern 


code 




KKb 




ObZ 










mesh 


32x16x32 


64x32x64 


128x64x128 


44x18x30* 










dz ■ 

'^''min 


0.03125 


0.015625 


0.0078125 


1.3 X 10-2 


17 


102 


3.2 X 10^ 


SL-3D 


Nu 


3.561 


3.631 


3.659 


3.61 


18 


32 


10^ 


SL-3D 


Nu 


3.119 


3.174 


3.197 


3.17 



Figure 2 



Fig. 2. Isothermal surfaces for several cases listed in Table 3 obtained by 
three-dimensional convection with strongly temperature-dependent viscosity sim- 
ilar to Ogawa et al. [34]. The calculations were carried out with 64 x 32 x 64 mesh 
divisions. 



mildly variable viscosity. 
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3.2.3 Benchmark for three-dimensional convection with strong viscosity vari- 
ation 

We also carried out the calculations similar to Ogawa et al. [34] in order to test 
our code for a much stronger temperature-dependence of viscosity. A steady- 
state convection in a box of 1.7 x 0.5 x 1 is considered. The impermeable 
and shcar-stress-free conditions are adopted along the all boundaries. We car- 
ried out calculations for several cases of Ogawa et al. [34] with several mesh 
divisions. In these calculations we varied Rat, the Rayleigh number defined 
with 77t, and the viscosity contrast = exp(£'T) between the top and bottom 
boundaries, as listed in Table 3. Figure 2 shows the convective flow patterns 
for several cases. We obtained the same flow patterns and the change in flow 
patterns depending on Rat and r,, as in Ogawa et al. [34]. The convective flow 
patterns are classifled into three-dimensional whole-layer convection (WL-3D) 
for Case 1 {Rat — 10^ and = 1), two-dimensional roll of whole- layer con- 
vection (WL-2D) for Case 4 {Rat = 10^ and = 10^), again in WL-3D for 
Case 16 {Rat = 10^ and T'f^ — 3.2 X lO"'), and in three-dimensional stagnant-lid 
convection (SL-3D) for Case 18 {Rat — 32 and — 10^). 

We summarize the comparison of the flow patterns and Nusselt numbers Nu 
obtained by the present code ("KKS") and those by Ogawa et al. [34] ("OSZ") 
in Table 3. The values of Nu obtained by both codes agree within at most 
4% deviation for all cases. The largest deviation for Case 1 may come from 
the differences in adopted mesh sizes. The present calculations are carried out 
with the minimum mesh size in z-direction of 1/128, which is more than four 
times flner than that employed in Ogawa et al. [34]. It is most likely that a 
flner spatial resolution is required for Case 1 than that employed in Ogawa et 
al. [34] in order to resolve the thermal boundary layers. 



3.3 Robustness and efficiency of multigrid iteration against spatial variation 
in viscosity 

To demonstrate the robustness and efficiency of the present algorithm against 
the spatial variation in viscosity, we performed convergence tests similar to 
Albers [22]. Here we solve Eqs. (1) and (2) only for prescribed distributions of 
buoyancy and viscosity. The distributions of buoyancy and viscosity are given 
by assuming the Rayleigh number Rat and the temperature-dependence of 
viscosity Et for a prescribed distribution of temperature. (The distributions of 
temperature used for the tests will be introduced below.) We take into account 
the change of the viscosity variation in the convecting vessel by increasing 
the global viscosity contrast = exp(i?T) from 1 to 10^°, and estimate a 
threshold value of (hereafter denoted by r^) below which the multigrid 
iteration converges. The values of velocity and pressure are initially set to zero. 
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Figure 3 



Fig. 3. Distributions of temperature used for convergence tests. Tliese temperature 
distributions were obtained by three-dimensional time-dependent simulation of ther- 
mal convection of fluids with temperature-dependent viscosity in a cubic box. The 
employed mesh division is 64 x 64 x 64. In (a) the several isothermal surfaces are 
shown, while in (b) the horizontally-averaged temperature (black lines) and the 
maximum of \grad{T)\ at height z (red lines) are plotted as a function of z. 

and are itcrativcly updated until the L^-norm of the residual r = b — Ax in 
Eq. (4) becomes smaller than that of b by eight orders of magnitude. All of 
the calculations were done with equally-spaced 64 x 64 x 64 mesh divisions 
and 6 grid levels, where the grid is coarsest for the grid level i — 1 and finest 
for i — 6. The mesh spacing is successively doubled as £ decreases, and the 
mesh spacing is 1/2 for grid level £ — 1. 

We used three distributions of temperature (hereafter denoted by Tempera- 
tures A, B, and C). In Figure 3 we present (a), in the left column, several 
isothermal surfaces and (b), in the right column, the plots of the horizontally- 
averaged temperature Th and the maximum of the magnitude of local temper- 
ature gradient \grad{T)\ma;^ at each height z, for these temperatures. The 
temperature fields are taken from snapshots of a three-dimensional time- 
dependent thermal convection of fluid with temperature-dependent viscosity 
in a cubic box. The Rayleigh number Rat defined by viscosity rjt is 10^, and 
Et = In(lO^) is assumed. In temperature A, which is taken from an early stage 
of the temporal evolution (nondimensional time t = 5.4 x 10~"^), the convective 
flow occurs only in the lowermost part of the box, and the local variation in 
the temperature is very small in the entire box (see the top panel of Figure 
3b). In temperature B (t = 1.52 x 10"^), the cold fluid with T < 0.5 sinks into 
the hot interior, and a large local variation in temperature occurs around the 
descending flow near the bottom surface (see the middle panel of Figure 3b). 
In temperature C {t — 3.6 x 10~^), the interior of the box is nearly isothermal 
and, as can be also seen from the smaller \grad{T)\jagj^ than in Temperature B, 
the temperature contrast between the descending flow and the surroundings 
becomes smaller. 

In addition, to flnd an appropriate implementation of the present algorithm 
for the multigrid operation, we employed four different types of smoothing 
procedures during the multigrid V-cycles. In these types, we varied the number 
of iterations A^^ by (6) for the pre- and post-smoothing calculations at each 
grid level, as listed in Table 4. In Type 0, Ns is taken to be 8 at all grid levels. 
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Tabic 4 

Numbers of pre- and post-smoothing calculations Ng on grid level i employed in the 
various types of the smoothing procedures in this study. Here A^grid is the number of 
employed grids and is the magnitude of global viscosity contrast. In the present 
calculations N„,.id is taken to be 6. 



Type 


Ns{£) on grid level £ 





8 


1 


8 X 2^sM-e 


2 


[ln(r,,) -Fl] X 8 X 2^sM-e 


3 


[In(r^) -h 1] X 8 X 4^g"d-^ 



Figure 4 



Fig. 4. Convergence behavior for various values of global viscosity contrast and 
for temperature fields in Figure 3 by using various types of smoothing procedures 
listed in Table 4. In (a) the number of V-cycles Ny is plotted against r^. In (b) 
the measure of total computational cost in the smoothing procedures Ws is plotted 
against r,,. The values are plotted only for < r^, below which the multigrid 
iteration converges. 

and is the smallest among the four types employed in these tests. In Type 1, 
Ng depends on the grid level i. The value of Ng is taken to be equal to that 
in Type at the finest grid = 6), and is successively doubled as the grid 
becomes coarser. In Types 2 and 3, in contrast, A^^, is assumed to be dependent 
on both the grid level £ and the global viscosity contrast r^. In both types, 
Ng at the finest grid is taken to be [In(r^) -|- 1] times larger than that in Type 
1, and the values of Ng is successively multiplied by two and four as the grid 
becomes coarser in Types 2 and 3, respectively. 



We compared the convergence behaviors of different types of smoothing proce- 
dures, for three distributions of T and various values of r^. In Figure 4 we show 
the plots against of (a) the number of V-cycles Ny and (b) the function Wg 
given by. 



6 



J2 n{i) X m{i) 

Ws = , (23) 

m[£ = 6) 



where n{£) is the total number of iterations by (6) in the smoothing proce- 
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Figure 5 



Fig. 5. Comparison in the evolution of the L^-norm of residuals on the finest grid 
(grid level £ = 6) during the initial five multigrid V-cycles between the case with 
the smoothing procedures of (a) Type 1 and (b) Type 0. Both calculations were 
performed with Temperature A and = 10*^. The red lines indicate the evolutions 
obtained by the smoothing calculations on the grid level £ = 6, while the blue arrows 
indicate the changes of residuals for ^ = 6 due to the coarse-grid correction. 

dures on grid level i, and m{i) = 4 x (64/2^^^)'^ is the number of elements of 
the unknown vector x (i.e., the number of unknown variables; see (6)) on grid 
level £. The numerator in (23) represents the total number of updating oper- 
ations of unknown variables during the entire multigrid iterations, while the 
denominator represents the number of updating operation during one smooth- 
ing calculation at the finest grid {i = 6). Namely, Ws is the measure of the 
computational cost of the multigrid iterations, and it indicates that the to- 
tal computational cost spent in the smoothing procedures during the entire 
multigrid iterations is equal to that virtually spent in l^g-times smoothing 
calculations on the finest grid. In the figure, the values are plotted only for 
T^T) 1^^, below which the multigrid iteration converges. 

The comparisons of for different types of smoothing procedures in Figure 
4 show that the value of becomes larger from Types to 3 for all tem- 
perature distributions. For Type 0, is 10^ for Temperature A, and 10^ for 
Temperatures B and C. By changing from Types to 1, increases to 10^ 
for Temperature A, and to 10^ for Temperatures B and C. The value of r"^ 
significantly increases by further changing to Types 2 and 3. For Type 2, is 
larger than 10^° for Temperature A, and 10^ for Temperatures B and C. For 
Type 3, rj^ is larger than 10^° for Temperature A, and 10^ for Temperatures 
B and C. From the comparison between these types, we conclude that the 
robustness of the multigrid iteration is improved by increasing A'^s- Iii partic- 
ular, a significant improvement is obtained for large by increasing Ng in 
proportion to [ln(rf,) -|- 1]. 

The multigrid iteration becomes more robust against for larger Ns, since 
larger number of iterations of (6) can reduce the errors more effectively. To 
see this more clearly, we show in Figure 5 the evolution of the L^-norm of 
residual during the multigrid iteration on the finest grid {£ — 6) for the cases 
(a) where the iteration converged with Type 1, and (b) where the iteration 
diverged with Type 0. For Type 1 where the multigrid iteration converges 
(Figure 5a), the L^-norm of residual on the finest grid is significantly reduced 
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after coming back from the smoothing on coarser grids. Even if the resid- 
ual is increased during the smoothing calculations on the finest grid (it can 
sometimes happen in this kind of decaying oscillation system; see eq. (15) 
or (19)), it is significantly reduced by the coarse-grid correction. As a result, 
the residual is successively reduced by the multigrid iterations. For Type 0, 
where the multigrid iteration diverges (Figure 5b), in contrast, the L^-norm of 
residual on the finest grid is increased after coming back from the smoothing 
on coarser grids. The smoothing calculations on the finest grid reduces the 
residual, except for the first multigrid iteration. However, the decrease in the 
residual by the smoothing calculations is too small to overcome the increase 
by the coarse-grid correction. As a result, the residual is successively increased 
by the multigrid iterations. From the comparison between these cases, we con- 
clude that sufficient number of smoothing calculations are necessary on coarser 
grids in order to ensure the convergence of multigrid iteration when a large 
spatial variation in viscosity is involved. 

Figure 4 also shows that the convergence rates of our method differ between 
the types of smoothing procedures. The comparison in the convergence rates 
between Types and 1 for small r,, indicates that a successive increase in 
A'^s to coarser grids significantly reduces the number of V-cycles Ny (Figure 
4a) as well as the computational cost Wg (Figure 4b). This stabihzing effect 
of increasing the number of smoothing iterations on coarser grids is in agree- 
ment with the features of previous results using multigrid methods [20,22, for 
example]. The comparison between Types 1 to 3 in Figure 4a indicate that 
Ny is significantly reduced by increasing Ns in proportion to [In(r^) + 1]. This 
feature is also consistent with the earher results [20,22, for example]. As can 
be seen in Figure 4b, however, the computational cost Wg of smoothing cal- 
culations increases as A^^ becomes larger from Types 1 to 3. This is because 
the cost of smoothing calculations during one V-cycle becomes significantly 
larger. 

We also note from Figure 4 that the convergence behavior differs between the 
prescribed temperature fields. This difference may reffect the difference in the 
magnitude of the local temperature gradient \grad{T)\ shown in Figure 3b. 
The different temperature distributions provide the different distributions of 
viscosity. The distribution of viscosity determines the nature of the coefficient 
matrix A in (4) and, hence, determines the convergence behavior of (6). In 
particular, a larger \grad{T) \ generates a larger local variation of viscosity for 
given r^, and makes the convergence of (6) more difficult. The convergence 
behaviors presented in Figure 4 are qualitatively consistent with the above 
conjecture. The multigrid iteration converges most easily for Temperature 
A with the smallest \grad{T)\, and least easily for Temperature B with the 
largest \grad{T)\. In Figure 4, the value of is the largest for Temperature 
A for any type of smoothing procedures. Although the values of are the 
same for Temperatures B and C when the same type of smoothing procedures 
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is employed, the number of V-cycles Ny is smaller for Temperature C than 
for Temperature B, except for the case with Type 2 for Temperature C and 
= 10®. These convergence behaviors are also consistent with those in the 
earlier studies [21,22, for example]. 



4 Discussion and concluding remarks 

We developed a numerical algorithm for solving mantle convection problems 
with strongly variable viscosity. Equations for conservation of mass and mo- 
mentum for highly viscous and incompressible fluids are solved iteratively by a 
multigrid method in combination with pseudo-compressibility and local time 
stepping techniques. In order to demonstrate its efficiency, the present algo- 
rithm has been implemented into a mantle convection simulation program 
based on the finite-volume discretization in a three-dimensional rectangular 
domain. Benchmark comparison with previous two- and three-dimensional 
calculations including the temperature- and/or depth-dependent viscosity re- 
vealed that accurate results are obtained even for the cases with viscosity 
variations of several orders of magnitude. Wc could also significantly improve 
the robustness of the numerical method against a spatial variation in viscosity 
by increasing the pre- and post-smoothing calculations in the multigrid oper- 
ations. We achieved convergence even for the viscosity contrasts up to 10^°, 
although the convergence rate deteriorates with increasing viscosity variations. 
The present algorithm can be further applied to the numerical models under 
more realistic conditions, such as spherical shell geometry and so on. 

The results of convergence tests described in Section 3.3 suggest that the 
convergence of multigrid method in combination with the present algorithm is 
determined by the accuracy of the coarse-grid correction when a large spatial 
variation in viscosity is incorporated. In the present convergence tests, the 
accuracy of coarse-grid correction was improved by increasing the amount of 
smoothing calculations on coarser grids during one multigrid iteration. This 
strategy always works, since the present iterative procedure never diverges as 
long as the local time stepping satisfies the CFL condition (see Section 2.3). 
However, it turned out that total computational costs significantly increased as 
the magnitude of viscosity contrast becomes larger (see Figure 4b) . Therefore, 
one needs to further accelerate the rate convergence of (6) in order to improve 
the efficiency of smoothing calculations and, in turn, to improve the accuracy 
of coarse-grid correction. 

The discussion in Section 2.3 also suggests that the rate of convergence of the 
present smoothing algorithm can be accelerated by choosing the matrix T in 
(6) more properly. As has been pointed out in Section 2.3, the matrix T acts as 
a "preconditioner" for (4). That is, the convergence of (6) becomes faster as T 
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better approximates A"^. In the present application to mantle convection 
problems, we had chosen T as a diagonal matrix by the use of the local 
time stepping approach. In addition, we had determined the values of nonzero 
elements of T from the assumption that both the rate of "diffusion" in velocity 
components and the rate of propagation of "pseudo-sound wave" were kept 
almost uniform in space. This choice of T corresponds to the preconditioning 
similar to the diagonal scaling of (4). However, apart from the local time 
stepping approach, we can use any arbitrary matrix T rather than a diagonal 
matrix. In other words, the convergence of (6) is accelerated if T preconditions 
(4) more properly than in the present paper. We can construct T, for example, 
by an incomplete factorization of A. Since the matrix T defined thus is most 
likely to be a better approximation of A~^ than that employed here, the 
convergence of (6) is expected to be faster than in the present cases. On 
the other hand, constructing T by an incomplete factorization has several 
disadvantages because (i) more memory is required to store all of the nonzero 
elements of T and (ii) the operation of an incomplete factorization is difficult to 
vectorize and parallelize. We should, therefore, take into account the specific 
computer architecture (such as scalar or vector processors) in choosing the 
most appropriate T in terms of preconditioning techniques. 

However, the most essential improvement possible for the present numerical 
method is the multigrid operation itself, rather than the smoothing calcula- 
tions during multigrid iteration. As has been suggested in Section 3.3, the 
reason why the multigrid iteration diverges for larger is that the coarse- 
grid correction becomes less efficient as increases (see Figures 4 and 5). 
The inefficiency of coarse-grid correction may come from an inappropriate 
treatment of the influences of variable viscosity in the multigrid operation. 
In the present numerical method, as well as in most of earlier models based 
on finite- volume discretization [26,20,22], we followed a "standard" strategy of 
multigrid. Namely, a linear interpolation was used for transferring the residual 
and error between adjacent grid levels. In addition, the discretized equations 
on coarser grids were derived by directly discretizing the differential equations 
on the particular coarser grids. However, as has been already acknowledged in 
the literature of multigrid [12,13], the standard strategy does not efficiently 
work when the differential equations contain strongly varying coefficients. One 
of the potential remedies is to develop the discretized equations based on 
Galerkin coarse-grid approximation [12,13], where the coefficient matrix A on 
a coarse grid is defined from the matrix on a fine grid A by, 

A = RAP, (24) 

where R and P are the restriction and prolongation operators, respectively. 
The Galerkin coarse-grid approximation can automatically define the coeffi- 
cient matrix A on the coarse grid, so as to sufficiently approximate A defined 
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on the fine grid. But it is not straightforward to use this approach in finite- 
volume models, since it makes the coefficient matrices on the coarse grids 
very complicated, with increased stencils, compared to those coming from the 
direct discretization. On the other hand, another potential remedy is to use 
the operator-dependent transfer operators [12,13]. These transfer operators 
utilize the structure of A with the variations in its coefficients and, hence, 
improve the accuracy of the restriction and prolongation operations. Indeed, 
Yang and Baumgardner [41] incorporated these transfer operators together 
with the Galerkin coarse-grid operators into their finite-element models, and 
demonstrated significant improvements in the robustness and efficiency of their 
multigrid procedures even for the cases with strongly viscosity variations. We 
expect that the robustness of our multigrid method presented in Section 3 
can be improved by incorporating the operator-dependent transfer operators 
into the present model, although their efficiency is left uncertain when used 
without the Galerkin approximation [22]. 

It is also an important issue to improve the robustness of the multigrid itera- 
tions against a "sharp" or "discontinuous" variation in viscosity, in addition to 
a "smooth" variation considered here. One of the major unsolved problems in 
the numerical study of mantle dynamics is to reproduce the motion of surface 
plates in the framework of mantle convection. Earlier numerical studies [42, 
for a review] had demonstrated that the generation of localized zones of low 
viscosity around the plate boundaries plays an important role in reconciling 
the fluid-like flow of mantle with the discrete motion of surface plates. How- 
ever, it is very difficult in essence for the multigrid method to deal with a local 
variation in viscosity as long as uniform mesh spacing is considered, since such 
local features are hardly "visible" on coarser grids. Recently, in several fields of 
geophysical fluid dynamics, new numerical techniques, such as the local mesh 
refinement [22], its combination with spectral-element method [43], and the 
wavelet-based method [44], are utihzed. Since these methods allow taking the 
spatial resolution non-uniformly in space, a local variation of viscosity can be 
handled accurately by, for example, using finer resolution around the regions 
of local variation. We thus speculate that the combination of these techniques 
together with the present method is one promising approach to reproduce the 
motion of surface plates in the numerical model of mantle convection. 
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Figure 1 



(a) Case la (b) Case 2 
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Figure 2 



(a) Case 1 {Rat = 10^ = 1) (b) Case 4 {Rat = 10^, = lO^) 
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Figure 3 



(a) isotherms (b) Th{z) and \grad(T)\^g,x{z) 

Temperature A {t = 5.4 x 10^'^) 
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Figure 4 



(a) number of V-cycle Ny (b) cost of smoothing procedures Ws 
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Figure 5 



(a) Type 1, Temperature A, = 10^ (b) Type 0, Temperature A, r.^ = 10^ 
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